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The dynamics of solitary gravity-capillary water waves propagating on the surface of 
a three-dimensional fluid domain is studied numerically. In order to accurately compute 
complex time dependent solutions, we simplify the full potential flow problem by taking 
a cubic truncation of the scaled Dirichlet-to-Neumann operator for the normal velocity 
on the free surface. This approximation agrees remarkably well with the full equations 
for the bifurcation curves, wave profiles and the dynamics of solitary waves for a two- 
dimensional fluid domain. Fully localised solitary waves are then computed in the three- 
dimensional problem and the stability and interaction of both line and localized solitary 
waves are investigated via numerical time integration of the equations. The solitary wave 
branches are indexed by their finite energy at small amplitude, and the dynamics of 
the solitary waves is complex involving nonlinear focussing of wave packets, quasi-elastic 
collisions, and the generation of propagating, spatially localised, time-periodic structures 
(breathers). 

1. Introduction 

Deep water gravity-capillary waves are relevant in a range of applications, including 
the understanding of the generation of waves by wind and the interpretation of satellite 
remote sensing data (see, for example I Zhang (1995) and Wright (1978) I. The study of 



solitary waves in this regime is of particular theoretical interest since it is only under the 
joint effects of surface tension and gravity that localised water waves in three-dimensions 
have been found in the water wave problem. Such localised waves have recently been 



observed in the experiments of Cho et al (2011) In the deep water limit, these Gravity- 



Capillary (GC) "wavepacket" solitary waves are of a fundamentally different nature 
from "long" solitary waves which are described approximately by the Korteweg-de Vries 
(KdV), Kadomtsev-Petviashvilli (KP) and related equations in shallow water. The later 
bifurcate from linear waves at zero wavenumber whereas the former bifurcate at a finite 
wavenumber, and hence their oscillatory nature. These oscillatory solitary waves were 



seen by Ablowitz & Segur (1979) and Akylas (1993) to be approximately described by 



solitary wave solutions of the focussing Nonlinear Schrodinger (NLS) equation which gov- 



erns the slowly varying envelope of monochromatic waves. In this regime, Akylas (1993) 
observed that if the phase speed (the speed of crests of the carrier wave) and the group 
speed (the speed of the envelope) are equal, the solitary waves of NLS will describe ap- 
proximately solitary waves in the primitive fluid equations. It is simple to show that at 
a minimum of the phase-speed, group and phase speed coincide, and that this occurs 
at nonzero wavenumber in GC waves for sufficiently deep water. (The condition is that 
the Bond number be less than 1 /3 which roughly corresponds to a depth greater than a 
centimetre in the air- water problem.) 
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In a two dimensional fluid domain (corresponding to a one dimensional free-surface and 
henceforth denoted as the 'ID problem'), wavepacket solitary waves were first computed 
in the full fluid equations by Longuet-Higgins (1989) and by Vanden-Broeck & Dias (1992) 



In three-dimensions (with a two-dimensional free-surf ace, denoted the '2D p roblem'), the 
first computations of steady solitary waves were by PAr^u et al (2005) with related 



work by Kim & Akylas (2005) and Milewski (2005) on reduced equations. Due to the 



highly oscillatory and spatially extended nature of these waves, accurate computations 
are challenging. Localised waves on a two-dimensional water surface are often called 
"lumps" , a name carried over from their shallow water counterparts which can be ap- 
proximated by the localised "lump" solutions of the KP equation, but we refer to them 
as wavepacket solitary waves. Q In this paper we focus only on the infinite depth case 
since for a water-air interface, any GC dynamics in water deeper than a few centimetres 
is essentially in the infinite depth regime. In this regime for an air-water interface, the 
waves bifurcate with a carrier wavelength of approximately 1.7 cm. and a speed of 23 
cm. /sec. 

The stability and dynamics of GC solitary waves in deep water has only recently been 



studied. In the ID problem, Calvo et al (2002) studied the question of linear stability and 
Milewski et al (2010) studied the time-dependent evolution of solitary waves (stability 
and collisions) in full potential flow. In 2D, there are far fewer studies. The transverse 
instability of line solitary waves (solitary waves of the ID problem trivially extended in 
the transverse variable) has been considered by Kim & Akylas (2007) and others. Line 



solitary waves are unstable to transverse perturbations of sufficiently long wavelength. 
Fully 2D dynamics have been considered by Akers & Milewski (2010) in a one-way sim- 



plified model and in a quadratic isotropic model in Akers & Milewski (2009) The main 



goal of this paper is to study 2D dynamics within a close approximation of the Euler 
equations. A simplification that we make, suggested and used in Kim et al (2012) for ID 
time-dependent solutions, is to take a cubic truncation of the scaled Dirichlet-to-Neuman 
operator that appears in the free-surface boundary conditions. (One important difference 
from the model of Kim et al (2012) is that we use the full surface tension term which 
results in a considerably better approximation of the full equations at larger amplitudes.) 
The approximation proposed is the simplest model that can quantitatively capture small 
and moderate amplitude Euler nonlinear CG solitary wave dynamics - a claim we support 
with ID comparisons. 

The focussing two-dimensional cubic NLS equation is central to the understanding 
of the existence and stability of these solitary waves. We shall see that this equation 
can be used to correctly predict the existence and certain instabilities of line and wave 
packet solitary waves, but does not capture the larger amplitude stability characteristics, 
the asymptotic dynamics of unstable waves nor the interaction of solitary waves. (The 
situation in ID is worse, since the NLS does not even capture the instabilities of arbitrarily 
small waves correctly - see for example, Calvo et al (2002) and Milewski et al (2010) ) 



This paper is structured as follows: in Section [2] we briefly present the derivation of 
the cubic truncation model we shall use and discuss what can be learned about the 
solitary waves from the associated NLS equation. In Section [3] we present the numerical 
results: ID comparisons between the cubic model and the full problem, followed by 
2D bifurcation diagrams, and stability and collision calculations. In Section U] we briefly 



f Naming these waves "wavepacket" solitary waves is more appropriate since it differentiates 
oscillatory solitary waves bifurcating from finite wavenumber from long waves bifurcating from 
zero wavenumber. 
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introduce variations on the cubic model, principally the addition of forcing and dissipation 
and treating the finite depth problem. 



2. Formulation 



2.1. Governing Equations 

Consider the three-dimensional free-surface water wave problem under the influence 
of both gravity and surface tension. Let (x, y) denote the horizontal plane, z the vertical 
direction and t time. The fluid is assumed to be inviscid and irrotational and therefore 
there exists a potential function 0, such that the fluid velocity (w, v, w) = {dxipi dy4'j dz4')- 
If the displacement of the water surface is designated by z = r]{x, y, t), then the governing 
equations for water waves read 



= 



Vt + iix4'x 



-> 
■Vy<l>y 



+ 



= 

+ = -V 
P 



z<7]{x,y,t) (2.1) 

as \x\ + \y\ -t- |z| ^ oo (2.2) 

at z — r](x,y,t) (2-3) 

at z = ?7(x,y,t) (2.4) 



where V = dy) is the horizontal gradient operator, and V- is the corresponding hori- 
zontal divergence operator. The constants g, p, a are the acceleration due to gravity, den- 



sity, and the coefficient of surface tension, respectively. Following Craig fc Sulem (1993) 
who worked in the canonical variables introduced by Zakharov (1968)[ the Kinematic 



and Dynamic boundary conditions (|2.3II2.4|) can be recast in terms of the free surface 
potential £^{x,y,t) — (j>{x^y,r]{x,y,t),t) and rj as 



Vt = G(77)e 

1 

6 



2(1-K iVryP) 

+(ve-v7?)2 



-r] + V ■ 



|VeP + 2(G(77)e)V^ Vt? - |VeP|V77|' 
V77 



(2.5) 



(2.6) 



G{ri) is a scaled Dirichlet to Neumann (DtN) operator yielding the vertical velocity of 
the free surface. It is defined by G{ri)^ — 4>z — 4>xrix — 'PyVy — vI+W'tF 'f'n, where 
(j)n is the derivative of the potential in the outward normal direction to the free surface 
and (j) satisfies (|2.m2.2p with (j){x,y,ri{x,y,t),t) — These equations have been also 

been nondimensionalized using a characteristic lengthscale ^ — [f^ ) i ^ timescale 

/ \l/4 / 2 

T — I \ , and a resulting velocity scale ^ = ( ^ ) ■ cgs units, g = 981cm/sec , 

and ^, the ratio of the surface tension coefficient and density, equals 73.5cm'^/sec^ for 
water. 



2.2. Expansion and Truncation 



Coifman & Meyer (1985) prove that if the L°°-norm and Lipschitz-norm of rj is smaller 



than a certain constant, then G is an analytic function of 77. It follows that the DtN 
operator can be naturally written in the form of Taylor expansion in rj, G = ^Gi. For 
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infinite depth, the first three terms of the Taylor series are given by 



Go - \D\ 

Gi{ri) = D- ^jD-GovGo 



G2{v) = -1[g, \D\' 



\D\^ if Go - 2Go ?7 Go 77 Go 



where D = — iV and \D\ = (— A)^/^. Using the cubic truncation of G and substituted 
into the kinematic and dynamic boundary conditions closes the evolution problem 



- m = (Gi 

- (1 - A)r/ = V • 

1 r 



G2)^ 



- V77 

2Go77Goe - 2rjA^) 



(2.7) 



(2.8) 



This formulation and approximation has reduced the three dimensional nonlinear wa- 
ter wave problem to a two dimensional one involving only the variables on the surface 
that is computationally, reasonably simple. In a doubly periodic setting, each term can 
be efficiently computed using a pseudospectral method and the fast Fourier transform 
(FFT). We shall henceforth call this model the cubic Dirichelet-to- Neumann Euler equa- 
tions (cDtNE) . Computational methods based on series truncations of the DtN operator 
have been used used in a variety of water wave problems and are summarised in de- 



tail by Nicholls (2007) For the ID problem where the full equations can be numerically 
integrated using a conformal map method we find (see the Results section) that the cu- 
bic truncation is extremely accurate. The next physically reasonable truncation for this 
problem, at fifth order, would involve substantially more computational resources. 

The cDtNE model can also be obtained from the fourth-order truncation of the kinetic 
energy part of the Hamiltonian expression of the surface water wave problem written in 
terms of the surface potential. The total energy of the fluid is the sum of kinetic and 
potential energies 



H =^ I dxdy 



[<Pl + 4>l + <i^'i)dz + ]^ j rfdxdy 



\J\ + |Vr/|2 - 1 ) dxdy. 



(2.9) 



and an approximate Hamiltonian can be derived by expanding the energy in powers of 
the ?7, ^. This takes the form 



1 

2 



e(Go + Gi + G2)i + -77' + ( + |Vr7p - "^dxdy (2.10) 



The equations (|2.7p and (|2.8I) can be expressed in canonical form in the sense of Zakharov (1968) 



r\t = ^ and — The system has further physical conserved quantities, of which 

mass and momentum 



Mass — 1] dxdy 



Momentum — / ^V77 dxdy 



(2.11) 
(2.12) 
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are used to monitor the global accuracy of numerical computations with our truncated 
equations. 

2.3. The Nonlinear Schrodinger Equation 

Traditionally, weakly nonlinear wavepackets are studied using the resulting cubic non- 
linear Schrodinger equation (NLS) for the modulational regime of monochromatic waves. 
It can be derived by substituting the ansatz: 

fA{X, Y, T)\ ^,(,.+,,_^,) ^ + ,2 j^^i V ^3 [^^2 V . . . (2.13) 



into (12. 5p . (|2.6p and ensuring that the series is well-ordered for t — 0{e^^). Here X = 
e{x — Cgt^, Y = ey and T = e^t where Cg is the group velocity in the wave propagating 
direction and "c.c." represents the complex conjugate of preceding terms. The wave 



envelope A can then be found to satisfy the NLS equation (Ablowitz & Segur (1979)) 

iAx + XiAxx + ^Ayy = fJ^\A\'^A (2.14) 

It is a trivial fact that substituting the same ansatz into the equations (|2.7p . (|2.8p instead 
yields an NLS equation with identical coefRcients - which is the motivation for choosing 
at least a cubic truncation of the DtN operator. We omit the details of the derivation, 
and just state the results. Choosing k = 1, Z = as the carrier wave, the phase and group 
velocity are equal, with the phase velocity at its minimum Cmin = V^. All of the waves 
we consider bifurcate from this point and exist only for c < Cmin ■ The coefficients of NLS 
are given by 

Ai = — , A. = — , A^ = -yV2. 

The solution to the original system is then as follows 

= eAe'^ - 2e^A'^e^'^ + ■■■ + c.c. (2.15) 

^ = -iV2eAe"^ + i2V2e^A^e^'^ + ■■■ + c.c. (2.16) 

where Q — x — \/2t. Since the NLS equation ()2.14p is of the elliptic or focussing type with 
solitary wave solutions in both one and two dimensions, and since the phase and group 
speed are equal at the chosen carrier wave, one can expect small amplitude solitary waves 
bifurcating from a uniform flow. These solitary waves of (j2.7l2.8p can be approximated 
by the NLS solitary waves found by solving the elliptic eigenvalue problem for p{x, y) 
and r2 

Ap + p^ = f7p, IIpII^^I, p-^Oatoo. (2.17) 

This is obtained by setting A = p{\\i\-'^/'^ X,\\2\-'^/'^Y)e-^^'^ . Denoting Ac = 

y/2 — c where c is the wave propagating speed, then, by the chosen scaling, Ac ~ ^le^ . 
The total energy of the Gravity- Capillary solitary wave bifurcating below the minimum 
phase speed is then calculated by (2.10). 

For the one dimensional problem (|2.17|) has the well known unique focussing NLS 
soliton, p — sech(x), — \, J = 2\/2 and thus 

||77|U«|-|'^'(Ac)'/' + |-|Ac (2.18) 
For the two dimensional case, which are the ones of interest here, there are countably 
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Figure 1. First three radially symmetric solutions of (|2.17[l . The corresponding energies of the 
wavepacket solitary waves for these three envelopes are 12.04, 79.39, 201.46 



many solutions to the problem (|2.17p (see Alfimov et al (1990) and references therein), 
and the first three radially symmetric solutions are shown in Figure [TJ For the ground 
state, called the Townes soliton in nonlinear optics (Chiao et al (1964) I, $7 w 0.204, 
J ~ 11.70, and, in the present hydrodynamic context. 



1/2, 



H 



4|AiA: 



11/2 



Ac 



P 



1/2 



12.04 



Ac 



(2.20) 



(2.21) 



The details of these computations can be found in Akers & Milewski (2010) An inter- 
esting particularity is that the physical energy for two dimensional wavepacket solitary 
waves is predicted to tend to a finite value of 12.04 as the amplitude approaches zero, 
which we shall verify in the cDtNE model. This is purely a consequence of scaling: the 
radially symmetric envelope's area is proportional to exactly countering the effect 
of decreasing amplitude of the wave whose energy density is e^. Shallow water lump 
solutions of KP (a valid approximation when the Bond number is greater than 1/3) do 
not have this property. The higher energy states of the eigenvalue problem ()2.17p can be 
associated with different families of travelling waves of the cDtNE model resulting in a 
remarkable "quantisation" of the energy of solitary GC waves in the 2D problem. 



3. Results 

The numerical solution of the cDtNE system is implemented on a periodic domain 
with Fourier pseudo-spectral methods, where all derivatives and Hilbert transforms are 
computed in Fourier space with spectral accuracy, while nonlinearities are computed 
pseudo-spectrally in real space. For traveling waves, the resulting algebraic system for 
the Fourier coefficients is solved using Newton's method using a monochromatic wave 
modulated by the solution to (|2.17p as initial guesses. The branches are computed 




Figure 2. Left: ID speed-amplitude bifurcation picture for elevation and depression solitary 
waves. The solid line, triangles, dotted lines and dashed lines correspond to cDtNE, Euler, 
leading order NLS and corrected NLS (from equation l2.19p respectively. Right: ID speed-energy 
bifurcation picture for elevation and depression solitary waves. The triangles denote the Euler 
solutions, the solid line is the cDtNE, and the dashed line is the prediction of NLS. 



through straightforward continuation methods. For time integration of the system, a 
classic fourth order Runge-Kutta method is used with the integrating factor method 



(see Milewski & Tabak (1999) for details) used to exactly integrate the linear part of 
the equation. The conserved quantities of the system are monitored and in all cases are 
preserved to a relative error of at most O{10~'^). All the computations are de-aliased 
with a doubling of Fourier modes. For two dimensional computations at least 256 x 64 
modes are used along the propagating and transverse directions respectively. It is often 
the small amplitude solutions that are most difficult to compute accurately since the 
spatial decay of those solutions is much slower. Thus, the computational the domain is 
gradually enlarged as the amplitude becomes smaller. 

3.1. The ID problem and model accuracy 
The bifurcation diagrams of ID Gravity-Capillary solitary waves for the full equations 



in deep water have been presented in Vanden-Broeck & Dias (1992) and others. In figure 



[5] we compare the speed-amplitude and speed-energy bifurcation diagrams of the cDtNE 
model to these and to the theoretical prediction of NLS. There, and elsewhere, we use 
either rj at the centre of the wave or the energy H as the amplitude parameter. The cDtNE 
model agrees well with the bifurcation picture for the Euler equations far beyond the NLS 
regime. Typical profiles of the depression and elevation branches of waves are show in 
figure [21 The model is remarkably accurate: at relatively large amplitudes which are far 
from the NLS regime the relative difference in profile between the full equations and the 
model are of order 10~^ and not visually discernible. Although a quantitative comparison 
of time-dependent dynamics is involved and beyond the scope of this paper (in particular 
methods used to solve the full equations usually do not use a uniform grid in a;), we show 
an example of the inelastic overtaking collision of two solitary waves computed with the 



full equations (reported in Milewski et al (2010) I and cDtNE truncation in FigurelH The 



results are extremely close and we see this as further evidence that the cDtNE model is 
an accurate representation of the Euler equations in a broad amplitude range of the CG 



regime. Comparisons that we have made with the cubic model of Kim et al (2012) show 



8 



Z. Wang and P. A. Milewski 





Figure 3. Typical free-surface solitary profiles of ID elevation and depression solitary waves 
for the cDtNE model far from the NLS wavepacket regime. Top: depression wave with 
?7(0) = -0.4406 and c = 1.3573. Bottom: elevation wave with ry(0) = 0.2055 and c = 1.3579. 
The solution to the full equations cannot be visually differentiated from these on this scale. For 
the depression wave, the maximum difference between the cDtNE and full solutions is 6 x 10~*. 




Figure 4. Overtaking collision of two depression waves of different amplitudes (minimum free 
surface heights of -0.36 and -0.12) shown in a frame of reference moving at the speed of the 
larger wave. From top to bottom: t = 0, t = 2000, t = 2500 and t — 3500. Only the larger wave 
survives the collisio n, with the smaller wa vepacket, to the left of it eventually dispersing. Left: 
Full equations from [Milewski et al (20To)] Left: cDtNE 



that retaining the full nonlinearity of the surface tension term is far more important than 
the nonlinearity associated with higher order corrections of the Dirichelet to Neumann 
map. One may conjecture that the surface tension term has a strong regularising effect 
that implies fast convergence of the DtN power series approximation throughout the 
evolution. 



3.2. The 2D bifurcation problem 

Two-dime nsional GC solitary w aves of the full equations in infinite depth were first 
computed in PArAu et al (2005) using finite differences and boundary integral methods. 




Figure 5. Left: Speed-amplitude bifurcation curves for elevation and depression solitary waves. 
The cDtNE, Euler, NLS and NLS plus correction are shown respectively by the solid, stars, dot 
and dash lines. Right: speed-energy bifurcation picture for both elevation (triangles pointing up) 
and depression (triangles pointing down). The bifurcation energy predicted by NLS is indicated 
by a square. 



A more resolved computation was performed by Parau and reported in Cho et at (2011) 
However, this later computation is, in our opinion, also under-resolved, particularly at 
small amplitudes where the waves are highly oscillatory and spatially extended, as it 
deviates considerably from both the cDtNE and the NLS results. Figure O shows the 
bifurcation diagram of GC solitary waves, as obtained from numerical solutions of the 



cDtNE, full potential flow (as computed by PSrSu and reported in Cho et at (2011) I, and 
both the leading order and leading order plus first correction of the NLS approximation. 
We conjecture that the cDtNE would be in quantitative agreement with full potential 
flow in the 2D problem given the accuracy of results at moderate amplitude in the ID 
case, and the agreement with NLS approximations at small amplitude. The discrepancy 
between the full Euler results and the expected amplitudes predicted by NLS had been 



noted in Cho et al (2011) Furthermore, we note that the cDtNE wave packets have finite 
energy at their bifurcation point, as predicted by NLS but which cannot be verified in the 
full Euler computations at current resolution (Parau, private communication). Typical 
profiles of the elevation and depression waves corresponding to the bifurcation curves 
presented in Figure [5] are shown in Figure [H 

All solutions presented so far are those whose profile envelope at bifurcation is de- 
scribed by the ground state eigenfunction of (|2.17p . In the 2D problem, we have also 
computed cDtNE solitary wave solutions whose envelope is approximated by radially 
symmetric higher modes. For example in Figure [12] we show the unstable evolution of 
the complicated travelling wave resulting from the solitary wave assigned to the second 
radial mode of (|2.17p (the wave is shown in the top- left panel of Figure [T^. This wave has 
an energy of 79.39 at bifurcation, much higher than the ground state. Non-radially sym- 
metric solutions to (|2.17|) are discussed in Alfimov et al (1990)[ but we did not attempt 



to compute solitary waves whose envelopes are not radially symmetric. 

3.3. Stability, focussing and wave collapse 

There are a few known results that guide us in a numerical study of 2D stability prob- 
lem. First, it is known that both depression and elevation plane solitary waves (ID waves 
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Figure 6. Top: representative of the depression wavepacket solitary wave solution to the 
cDtNE equation with r;(0, 0) = —0.41 and c — 1.3907. Bottom: representative of the elevation 
wavepacket solitary wave solution to the cDtNE equation with 77(0, 0) = 0.23 and c = 1.3842. 



extended in the second dimension) are linearly unstable with respect to sufficiently long 
small perturbations in the transverse direction. This has been shown both within an NLS 



approximation Rypdal fc Rasmussen (1988), and using arguments based on Unearization 
of the full equations Kim & Akylas (2007) Second, in 2D the underlying focussing NLS 



equation (|2.14p is well known to exhibit a finite-time focussing blowup called wave col- 
lapse when, in an unbounded setting, the initial conserved energy E = J | VAp — ^ |^|^, is 



negative (see Zakharov (1972) and Sulem & Sulem (1999) I. (Note that this energy is not 



the same one arising from the cDtNE equations.) The result is obtained by a virial argu- 
ment, whence Mq = J |Ap is conserved and M2(i) = J{X^ +Y'^)\A\'^ satisfies M" = 8E. 
Solitary waves which are solutions to the eigenvalue problem ()2.17|) are critical with 
E = 0. Thus, a small negative energy perturbation leads to focussing and blowup and a 
positive energy perturbation leads to spreading of the underlying wave envelope. We thus 
expect, and will confirm numerically, the instability of localised solitary wave solutions 
for nonzero envelope energy perturbations in the near-NLS limit of the cDtNE model 
(where solitary wave envelopes are well approximated by the NLS). We note that the 
simple virial argument is not available in a periodic setting (see Sulem & Sulem (1999) I 
but nevertheless seems to predict stability very well in our computations on periodic 
domains. Third, both elevation and depression wave branches have critical point in their 
speed-energy relation (see Figure [5]), which can lead to an exchange of linear stability 
of the eigenfunction of the linearized problem associated to the translational invariance 



symmetry (see, for example, Akylas & Cho (2008) I. This is the instability commonly 



observed in small amplitude elevation waves of the ID problem leading to the eventual 



development of a depression wave Milewski et al (2010) 



All figures of solitary wave dynamics presented are shown is a frame moving with the 
speed of the wave that was used construct the initial data. In this frame the main features 
of the wave evolve slowly. 



For plane solitary waves, the linear analysis based on NLS (see Akers & Milewski (2010) ) 
shows that the transverse perturbation e*'^ is unstable, when the wave number I in the 
y direction satisfies 



1^1 < 



3/i 



8A2 



1/2, 



(18)i/4(Ac)i/2 A 



(3.1) 



We confirm this in our numerical experiments where a plane solitary wave is perturbed 



Gravity- Capillary Solitary Waves 11 




100 200 300 400 500 600 700 800 

Time 



Figure 7. Transversal instability of the plane solitary wave, c = 1.3994, r;(0) — —0.247. 
The transverse perturbation is obtained by taking an initial data of the form 
rj{x,y,0) = [1+0.1 cos(aZcJ/)]'7(a;), where 7)(a;) is a travelling wave solution in ID, Ic = 18^/* Ac^/^, 
and a = 0.8 (solid line), a = 0.9 (dashed line), a = 1.1 (dotted line), a = 1.2 (plus line) 
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Figure 8. Evolution of the transverse instability to a localised solution in the case of a = 0.8 
in Figure [7] The contours are shown from top to bottom at t = 0, t = 125, t — 175. 



with four transverse perturbations of different wavelengths and the subsequent nonhn- 
ear evolution is compared. The perturbation wavenumbers are / = 1.2lc, l.l/c,0.9Zc, and 
0.8/0 respectively. As shown in Figure [7] the first two perturbations do not destabilise 
the plane wave whereas the next two do. The subsequent evolution of the instability 
shows a focussing behaviour reminiscent of the underlying collapse dynamics of NLS as 
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shown in Figure HI This intermediate time focussing is arrested by the generation of a 
travelhng breather: a propagating periodic-in-time focahsed structure. This structure is 
best described as a locahsed depression sohtary wave with periodic amphtude modula- 
tion, and which appears to be stable (see evolution for i > 175 in Figure [7]) in all our 
experiments. These breathers are very common and also occur in the nonlinear evolution 
of instabilities of small amplitude fully localised solitary waves below. 

Throughout our computations in this section, instabilities and wave interactions will 
invariably lead to some high frequency dispersive radiation and thus, due to our use of 
a periodic domain, the remaining coherent structures are embedded in a "sea" of linear 
ripples. This is visible in most computations and can also be seen as further evidence of 
the stability of the resulting structures. 

The stability of localised traveling waves is considered next. First, the virial argument 
sketched above implies that near the bifurcation point, both elevation and depression 
solitary waves are linearly unstable. The virial argument would predict eventual blowup 
for negative energy perturbations and this is not observed, in all cases the blowup is 
arrested by the generation of a larger amplitude breather. Dispersive spreading consistent 
with the envelope spreading predicted by the virial argument is observed for positive 
energy perturbations. In the computations that we present we perturb the exact solitary 
wave solution by a small multiple of itself, by taking initial data r]{x,y,0) = (1 + S)rj 
where rj is the computed travelling wave. Since, for the perturbed wave, 



E = J \VA\'-^\A\*^-5 1 \A\' 



where A is the envelope of the solitary wave, negative energy perturbations correspond 
to 5 > and positive energy ones to S < 0. In the left panel of Figure IH] we show the 
typical evolution of the amplitude for a perturbed small amplitude depression solitary 
wave. For a negative energy perturbation we see focussing arrested by the formation of 
a breather whereas for positive energy the amplitude decreases monotonically as a result 
of dispersive spreading. The case shown is for a depression wave, however the elevation 
wave dynamics is broadly similar. In Figure [TU] we show four snapshots of the breather 
evolution resulting from the unstable small amplitude depression wave at later times. 

For a large class of problems, linear stability may change at a critical point of H{c), the 
speed-energy curve. This necessary condition for instability was observed by Saffman (1985)[ 



for gravity waves and has since been extended to many other situations. Akylas & Cho (2008) 
present an application of this result to a wavepacket solitary wave in a model equation. 
In the present problem, both depression and elevation waves have critical points in H{c). 
We denote solitary waves with speed lower than this critical speed "large amplitude" and 
those with speed larger than it "small amplitude" since in all our computations amplitude 
is a monotonic function of speed. Small amplitude waves' bifurcation diagram and linear 
stability are well described by the NLS equation whereas large amplitude waves' are not. 
Our numerical simulations show that for depression waves, the exchange of stability does 
take place at the minimum point in H[c). The evolution of large amplitude depression 
waves' amplitudes when subject to perturbations is shown in the right panel of Figure 
[9l The waves are stable regardless of the sign of 5. Solitary waves (or breathers which 
are small perturbations of the solitary waves) propagate in the midst of small linear 
dispersive waves that have been shed by the perturbed initial data. 

All our computations show that elevation waves remain unstable at large amplitude. 
In Figure [Tl] an unstable large elevation solitary wave subject to a smallnegative energy 
perturbation and evolves into a depression breather. Similarly, positive energy perturba- 
tions will also yield depression breathers (whose energy is much smaller than the elevation 
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Figure 9. Left: evolution of maximum trough depth for perturbed depression waves in the small 
amplitude regime (77(0,0) = —0.24 and c — 1.4103). Negative energy perturbation {S = 0.01) 
is shown by a solid line, and positive energy {S = —0.01) by a dashed line. Right: evolution 
of maximum trough depth (normalised by initial maximum trough depth) for depression waves 
in the large amplitude regime (r;(0, 0) — —0.49 and c — 1.3873). Positive energy perturbation 
{5 = 0.01) is shown by a solid line, and negative energy [S = —0.01) by a dashed line. 




Figure 10. From top-left to bottom-right: snapshots of one period of a breather taken at 
t = 1240, 1440, 1520, 1620 for the small amplitude negative energy computation discussed in 
Figure ^ 
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Figure 11. From top-left to bottom-right: snapshots of the evolution of a large amphtude 
unstable elevation wave with ??(0,0) = 0.2196, c = 1.3907 in cDtNE at t = 0,120,200,270. 
The perturbation is 0.01 times the original solitary wave. Some waves are radiated during the 
transition and the unstable wave evolves into a breather. 



waves (see Figure [5]). The instability is initially similar to that of a ID depression wave 



(see Milewski et al (2010)): a symmetry breaking whereby the leading trough grows at 
the expense of the trailing one (this instability is associated to the translational invariance 
mentioned above). What follows is collapse focussing that is arrested by the formation 
of a depression breather. 

From our calculations, we believe there are families of periodic (breather) solutions of 
different periods and amplitude for each fixed energy above the minimum of the depres- 
sion solitary waves. The orbits of these breathers in phase space are centred around the 
fixed point of stable depression solitary waves and the precise bifurcation diagram for 
them would require computing exact periodic localised structures which is beyond the 
scope of this paper. A summary of the stability results in this section is shown in table 



Lastly we compute the evolution of a complex solitary wave of the cDtNE equations 
corresponding to a higher mode of the nonlinear eigenvalue problem discussed previously, 
and this is shown in Figure fT2l In this case all perturbations triggered rapid instabilities, 
but the nonlinear evolution shows a remarkable dynamics with eventual focussing into 
four depression breathers of different amplitudes. 

3.4. Collisions 

Given the stability of large amplitude depression waves, we have numerically computed 
their collisions. We have only computed head-on and overtaking collisions of pairs of 
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positive energy perturbation negative energy perturbation 

small depression unstable— >-large breather unstable^disperses out 

small elevation unstable^large breather unstable^disperses out 
large depression stable stable 

large elevation unstable— ^large breather unstable^large breather 

Table 1. Summary of stability Results 





-60 -40 -20 20 40 60 -60 -40 -20 20 40 60 

Figure 12. Instability of a higher energy solitary wave with ri{0,0) = —0.4028, c = 1.4061. 
Snapshots of the evolution are shown at t = (top left), t — 800 (top right), t — 1200 (bottom 
left), t = 1600 (bottom right). Four depression breathers emerge from the complex behaviour. 



waves although in a 2D problem there is a wide range of possible collision scenarios. 
For head-on collisions, the interaction time between the two waves is insufhcient for any 
strong nonlinear effect to take place and we have only observed very small oscillations 
that result from the small inelasticity of the collision. The waves essentially traverse each 
other. The more interesting case is the overtaking collision. Here, the small difference in 
solitary wave speeds implies that the calculations must be carried out over long times. In 



the ID problem collisions were of two types of inelastic collisions Milewski et al (2010) 



collisions where both waves survived and collisions where only the larger wave survived 
when their amplitude difference was large. Here, we have only observed quasi-elastic 
collisions where both waves survive and the primary effect of the overtaking collision is a 
rapid and large phase shift of the order of one envelope wavelength. Figure [13] shows the 
before and after free surface profiles of the waves and Figure [14] shows the resulting waves' 
trajectories in {x,t) space. The collisions are weakly inelastic and the wave amplitude is 
mildly attenuated with some dispersive radiation present during the collision. 
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Figure 13. Collision of two stable depression solitary waves travelling in the positive x-di- 
rection. Left: Initial data consisting of a shifted superposition of two waves with with 
?7(0,0) = -0.56, c = 1.3782 and 77(0,0) = -0.49, c = 1.3873. Right: Solution after interaction 
at t = 4000. Note that the collision is not completely elastic, with their amplitudes decreasing 
slightly as a result of the collision. The solution was computed in a frame of reference moving 
at speed 1.3828. 
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Figure 14. Filled contour plot of the centreline of the dynamics ri{x,0,t) for the collision 
reported in Figure [131 Note that the waves' core remain distant during the collision and that 
after the collision the waves have exchanged identities and undergone a phase shift. 



4. Extensions 

There are several extensions of the model which, while the detailed study is beyond the 
scope of this paper, may be useful in studying problems of this type. We have considered 
fluids of infinite depth only (which in the case of GC flows in water is a good approxi- 
mation for depths exceeding a few centimetres). If the effect of finite depth is required, 
the formulation can be modified simply by taking Go(^) = 1^1 tanh(/i|Z?|) where h is the 
mean depth of the fluid normahsed by the capillary gravity length scale and is inversely 
related to the Bond number. In particular the shallow water lump dynamics modelled by 
the Kadomtsev-Petviashvilli (KP) equation should be recovered when the Bond number 
is greater than 1/3. 

Given the small length scales of the waves, it may also be of practical interest to add 



the effects of viscosity and forcing. In the analysis of their experiments Gho et al (2011) 
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adjusted a simple one-way model proposed in Akers & Milewski (2009) by adding forcing 
and a viscous damping. The more accurate cDtNE equations can be modified to include 



small viscous damping effects by using the approximation presented in Dias et al (2008) 



whereby dissipation is modelled through the modification of both kinematic and dynamic 
boundary conditions: 



Vt 

6 + (1 ■ 



A)r, 



2i?e-iAr;+(Gi+G2)e 
2Re-'^A(, + P{x,y,t) + \7 

1 



(GoO(Go^ - 2Go?7Go^ - 2r;AC) - |VC| 



V?7 



(4.1) 



(4.2) 



Here, P{x, y, t) is the pressure forcing and the Reynolds number, which controls the 
dissipation rate, is given by 

V 



Re = 



where v is the kinematic viscosity of the fluid. For the case of an air- water interface in the 
regime that we considered in this paper, we obtain a Reynolds number of approximately 
500. 

For small amplitude waves, the variation in the transverse direction is not significant 
compared to that in the propagation direction. One can therefore propose to assume 
this a priori and simplify the system by deleting all the nonlinear terms which include y- 
derivatives obtaining a "weakly transversal" model which still has the same NLS equation 
describing wave packets as the full problem. The Hamiltonian for this model reads 



H[ri.£]^ I ^e(Go + Gi + G2)e 



1 



1 



(Vl + ^x) dxdy 



(4.3) 



where 



Go — ( — dxxj 

Gi = -dxrjdx 
_ 1 

G2 



1/2 



^Gorfdxx + ^dxxV^Go + Gor]GQT]Go 



(4.4) 



Calculations performed with this approximation show excellent agreement with cDtNE 
at small amplitudes and only qualitative agreement for larger amplitude solitary waves. 



5. Conclusions 

The dynamics of Gravity-Capillary solitary waves on the surface of three-dimensional 
fluid was studied. The only approximation made was a cubic truncation of a scaled 
Dirichelet to Neumann map that provides the normal velocity of the free-surface given 
its tangential velocity. In two dimensions, where comparisons to the untruncated problem 
(i.e. fully nonlinear free-surface potential flow) can be accurately measured, this trun- 
cation is remarkably accurate in modelling small and moderate amplitude waves. We 
conjecture that the same is true in three-dimensions. 

There are undoubtedly in this problem infinitely many branches of solitary travelling 
wave solutions as this is the prediction of the associated envelope NLS analysis about 
the bifurcation point. We compute examples of three of them: elevation and depression 
solitary waves arising from the simplest eigenfunction of the NLS problem and a more 
complex wave arising from a higher eigenfunction of the same problem. The localised 
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solitary waves have the surprising property that as their amphtude (from peak to trough) 
decreases to zero as the bifurcation point is approached, their physical energy tends to a 
finite positive value, quantised by the different eigcnfunctions of the NLS equation. 

The instability and subsequent evolution for one dimensional line solitary waves and the 
various two dimensional solitary waves have been explored numerically by perturbing the 
waves and computing the solution through accurate pseudo spectral based methods. All 
solitary waves are found to bc^ unstable with the notable exception of larger amplitude 
depression waves. These waves together with travelling breathers, which are periodic 
travelling cycles oscillating about these travelling states, are stable and appear to be 
attractors in the long time evolution of the problem. 

The focussing NLS equation adequately predicts the bifurcation and linear stability 
properties of small amplitude solitary CG waves. The collapse singularity of initial data, 
however, is not observed in our computations leading to the conjecture that an appropri- 
ate envelope model for this problem is a cubic-quintic NLS equation where the quintic 
term is defocussing. 
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